Linear response theory for long-range interacting systems in quasistationary states 
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Long-range interacting systems, while relaxing to equilibrium, often get trapped in long-lived qua- 
sistationary states which have lifetimes that diverge with the system size. In this work, we address 
the question of how a long-range system in a quasistationary state (QSS) responds to an external per- 
turbation. We consider a long-range system that evolves under deterministic Hamiltonian dynamics. 
The perturbation is taken to couple to the canonical coordinates of the individual constituents. Our 
study is based on analyzing the Vlasov equation for the single-particle phase-space distribution. The 
QSS represents a stable stationary solution of the Vlasov equation in the absence of the external 
perturbation. In the presence of small perturbation, we linearize the perturbed Vlasov equation 
about the QSS to obtain a formal expression for the response observed in a single-particle dynam- 
ical quantity. For a QSS that is homogeneous in the coordinate, we obtain an explicit formula for 
the response. We apply our analysis to a paradigmatic model, the Hamiltonian mean-field model, 
which involves particles moving on a circle under Hamiltonian dynamics. Our prediction for the 
response of three representative QSSs in this model (the water-bag QSS, the Fermi-Dirac QSS, and 
the Gaussian QSS) is found to be in good agreement with A''-particle simulations for large A^. We 
also show the long-time relaxation of the water-bag QSS to the Boltzmann-Gibbs equilibrium state. 



PACS numbers: 61.20.Lc, 05.20. Dd, 64.60. De 

I. INTRODUCTION 

Systems with long-range interactions are ubiquitous 
in nature [TH3]. In these systems, the interaction po- 
tential between two particles decays asymptotically with 
separation r as 1/r", where a is smaller than or equal 
to the spatial dimension of the system. Examples are 
dipolar ferroelectrics and ferromagnets, self-gravitating 
systems, non-neutral plasmas, two-dimensional geophys- 
ical vortices, etc. Long-range interactions result in non- 
additivity, thereby giving rise to equilibrium properties 
which are unusual for short-range systems, e.g., a nega- 
tive microcanonical specific heat, inequivalence of statis- 
tical ensembles, and others [T]. 

Long-range systems often exhibit intriguingly slow re- 
laxation toward equilibrium [3HZ]- Such slow relaxation 
has been widely investigated in a model of globally cou- 
pled particles moving on a unit circle and evolving un- 
der Hamiltonian dynamics. This model is known as the 
Hamiltonian mean- field (HMF) model [5]. In this model, 
a wide class of initial conditions relaxes to equilibrium 
over times that diverge with the system size. It has 
been demonstrated that, depending on the initial con- 
dition, the relaxation to equilibrium for some energy in- 
terval occurs through intermediate long-lived quasista- 
tionary states (QSSs). These non-Boltzmann states in- 
volve slow evolution of thermodynamic observables over 
time, and have a lifetime which grows algebraically with 
the system size An immediate consequence is that 
the system, in the limit of infinite size, never attains the 
Boltzmann-Gibbs equilibrium, but remains trapped in 
the QSSs. Generalizations of the HMF model to include 



anisotropy terms in the energy , and to particles which 
are confined to move on a spherical surface rather than 
on a circle [TUJ have also shown slow relaxation toward 
equilibrium and the presence of QSSs. 

Dynamics of systems with long-range interaction, in 
the infinite-size limit, is described by the Vlasov equa- 
tion that governs the time evolution of the single-particle 
phase-space distribution [TT]. This equation allows a 
wide class of stationary solutions. Their stability can be 
determined from the temporal behavior of small fluctua- 
tions by linearizing the Vlasov equation around the sta- 
tionary solutions. Stable stationary states of the Vlasov 
equation correspond to QSSs of the finite-size system. 
The first clear demonstration of this correspondence was 
achieved for the HMF model [1]. In this paper, we refer 
interchangeably to stable stationary states of the Vlasov 
equation and QSSs of the finite-size dynamics. 

The study of the time evolution of fluctuations around 
a stationary solution of the Vlasov dynamics is relevant 
in many applications, for instance, in the phenomenon of 
Landau damping that arises due to energy exchange be- 
tween particles and waves in an electrostatic plasma [H] . 
In this process and in many others, spontaneous statis- 
tical fluctuations around a stable stationary state of the 
Vlasov equation are considered, and their rate of expo- 
nential decay in time is determined by solving an initial 
value problem involving the linearized Vlasov equation 
through a careful use of the Laplace-Fourier transform. 

In this paper, we follow a different approach to the 
study of fluctuations by the linearized Vlasov equation. 
Inspired by Kubo linear response theory [13] , we analyze 
the response of a Vlasov-stable stationary state to the 
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application of a small external perturbation described by 
a time-dependent term in the Hamiltonian. The pertur- 
bation induces forced fluctuations around the stationary 
state that we treat to linear order in the strength of the 
perturbation, and study their evolution in time by using 
the linearized Vlasov equation. Such forced fluctuations 
are known to be generically flnite for Boltzmann-Gibbs 
equilibrium states. We show here theoretically that they 
are finite and small, of the order of the perturbation, also 
for Vlasov-stable stationary states. We support our anal- 
ysis with iV-particle numerical simulations of the HMF 
model. 

The paper is organized as follows. In Sec. we de- 
velop the linear response theory for a general QSS by 
using the Vlasov framework. In Sec. |III[ we specialize 
to a QSS that is homogeneous in the coordinate, and de- 
rive a closed form expression for the change induced by 
the external perturbation in a single-particle dynamical 
quantity. Section |IV] is devoted to the application of the 
theory to study the response of three representative ho- 
mogeneous QSSs in the HMF model, namely, the widely 
studied water-bag QSS, the Fermi-Dirac QSS and the 
homogeneous equilibrium state, which is also a QSS. In 
the following section, we compare results from A'^-particle 
numerical simulations of the HMF dynamics with those 
from the linear response theory, and obtain good agree- 
ment. We also discuss the long-time relaxation of the 
water-bag QSS to Boltzmann-Gibbs equilibrium under 
the action of the perturbation. We draw our conclusions 
in Sec. ED 

II. LINEAR RESPONSE THEORY FOR QSS 

Consider a system of N particles interacting through 
a long-ranged pair potential. The Hamiltonian of the 
system is 

N 2 1 ^ 

where qi and pi are, respectively, the canonical coordi- 
nate and momentum of the ith particle, while v{qi — qj) 
is the interaction potential between the ith and jth parti- 
cles. We assume that v{q) is suitably regularized for zero 
argument. We regard q^'s and pi's as one-dimensional 
variables, though our formalism may be easily extended 
to higher dimensions. The mass of the particles is taken 
to be unity. The factor 1 /N in Eq. ([T]) makes the energy 
extensive, in accordance with the Kac prescription |14| . 
In this work, we take unity for the Boltzmann constant. 
The system evolves under deterministic Hamiltonian dy- 
namics: the equations of motion for the ith particle are 

d 1 ^ 

q^=P^, P^ = -Q^ -^^v{q,-qJ), (2) 
i<j 

where dots denote differentiation with respect to time. 



We start with the system in a quasistationary state 
(QSS) at time t = 0, and apply an external field K{t). 
A QSS represents a stable stationary solution of the dy- 
namics ([2]) in the limit N 00. For finite N, however, 
size effects lead to instability and a slow relaxation of 
the QSS to the Boltzmann-Gibbs equilibrium state over 
a timescale that diverges with TV [WTligifTU]. 

Assuming the field K{t) to couple to the coordinates 
of the individual particles, the perturbed Hamiltonian is 

N 

Hit) ^Ho + i/ext =Ho- K{t) h{q^). (3) 

1=1 

Here, b{qi) denotes the dynamical quantity for the zth 
particle that is conjugate to K{t). The equations of mo- 
tion are modified from Eq. ([2| to 

q^=P^, P^ = -|:^E«fe-9.) + ^W^- (4) 

In this work, we study the temporal response of the 
initial QSS to the field K{t), in particular, the linear 
response. We ask: How does a single-particle dynamical 
quantity a{q), that starts from a value corresponding to 
the QSS, evolve in time under the action of K{t)l We 
seek answers to this question by considering the system 
in the limit — > 00, so that size effects are negligible and 
the evolution of the QSS is due to the field K{t) alone. 
We also regard K{t) to satisfy the following conditions: 
K{t) is a monotonically increasing function of t and has a 
value <C 1 at all times, K{t = 0) = 0, and K{t — 00) = 
a constant much smaller than 1. While discussing the 
time-asymptotic response, we will mean the ordering of 
limits A^ — 00 first, followed by t — 00. Note that 
the perturbed dynamics Eq. Q does not conserve the 
total energy of the system as does Eq. ([2| , although the 
variation is expected to be small for small K(t). 

The framework we adopt to address our queries is that 
of the Vlasov equation for the time evolution of the single- 
particle phase space distribution. For a system like Eq. 
([1]) in the limit N 00, such an equation faithfully 
describes the A^-particle dynamics in Eq. ([2| [TTJ \T^. 
That for small K{t) the Vlasov equation describes the 
perturbed dynamics of Eq. Q in the infinite-size limit 
is illustrated later in the paper by comparing the predic- 
tions of our analysis with A^-particle simulations for large 
N. 

Let the function fd{q,p,t) count the fraction of parti- 
cles with coordinate q and momentum p at time t: 

1 ^ 

Uq,P,t)^-J2s{q- q, {t))s (p - p, (t)) . (5) 

In the limit A^ — >■ 00, the function fdiq,P,t) converges to 
the smooth function f{q,p,t), which may be interpreted 
as the single-particle phase space distribution function. 
The Vlasov equation for the time evolution of f{q,p,t) 
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may be derived by using the equations of motion Q and 
following standard approaches jlH I15j. and is given by 



dl 
dt 



-C{q,p,t)[f]f = 0, 



(6) 



where the operator C{q,p, t)[f] is given by 



while $((7,t)[/] is the mean-field potential: 

'^i<l,tm^ JJ dq'dp' v{q-q')f{q',p',t). (8) 

We investigate the response of the system to the ex- 
ternal field by monitoring the observable 



a('7)>(^)= // dqdp a{q)f{q,p,t). 



(9) 



To obtain its time dependence, we need to solve Eq. ([6| 
for f{q,p,t), with the initial condition 



(10) 



Here, fa{q,p) characterizes a QSS, i.e., a stable station- 
ary solution of the Vlasov equation for the unperturbed 
dynamics Thus, fo{q,p) satisfies 



where 



and 



'^^o(g,p)[/o] = -P 



>Co(g,p)[/o]/o-0, 

d , 9$(g)[/o] d 



dq dp' 



dq'dp' v{q-q')h{q',p'\ 



(11) 



(12) 



(13) 



To solve Eq. ([6j) for K{t) <C 1, we expand f{q,p,t) to 
linear order in Kit) as 

/(g,P,t) = /o(9,p) + A/(g,p,i), (14) 

with the initial condition 

A/(g,p,0) =0. (15) 



Substituting Eq. (14 1 in Eq. and separating terms 
to order 1 and K{t), we get, respectively. 



dfo 
dt 



>Co(g,p)[/o]/o = 0, 



(16) 



and 



dAf 
dt 



- Co{q,p)[fo]Af = C,^t{q,P,t)[Af]fo. (17) 



Here, the operator 

. , Nr. .1 a$(q,t)[A/] d ^dh d , 

/:e.(.,P,0[A/] = ^^^-i^(t)^^ (18) 



describes the effects of the external field, which are two- 
fold: (i) to generate a potential due to its direct coupling 
with theparticles, and (ii) to modify the mean-field po- 
tential (|8| from its value $(g)[/o] in the absence of the 
field. Defining an effective singlc-particlc potential, 

Ueff((?,i)[A/] =$(q,t)[A/] -X(t)6(g), (19) 



Eq. (17) may be written as 



9A/ 9«eff(g,i)[A/]a/o 
i'o(g,p)[/oJA/ = - 



dt 



dq 



dp 



(20) 



Equation ( 16 ) is satisfied by virtue of the definition of 



foiq,p)- We thus solve Eq. (20) for Af{q,p,t) in order 
to determine f{q,p,t) from Eq. (14). With the condition 
(15), the formal solution is 



Af{q,p,t) 



Z 

J dT e(*- 



-r)£o(<7,p)[/o] ^^cff(g,T)[A/] dfo{q,p) 



dq dp 


(21) 

Using Eq. (21) in Eqs. ([9| and (14) gives the change 
in the value oi(a{q)){t) due to the external field: 



{Aaiq)){t)= II dqdpa{q)(^fiq,p,t)- fo{q,p) 



= / dr 



dqdp a((j)e(*-^)^''(«'P)[^ol 



dv,ft{q,T)[Af] dh(q,p) 



dr 



dp 

da{t~T) dv,s{q,T)[Af]\ 



dp 
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/o 



(22) 



Here, angular brackets with /g in the subscript denote 
averaging with respect to fo{q,p), e.g.. 



while 



a(g))/o = / / dqdp a{q)foiq,p), 



a{t-T) = e-(*-^)^°(''f)[^''la(9) 



(23) 



(24) 



is the time-evolved a{q) under the dynamics of the unper- 
turbed system. In obtaining the last equality in Eq. (22 ), 



we have used the definition of £oj have performed integra- 
tion with respect to g, and have assumed the boundary 
terms involving fo{q,p) to vanish. 

Defining the Poisson bracket between two dynamical 
variables g{q,p) and g'{q,p) in the single-particle phase 
space as 



{9{Q,p),9'{q,p)} 



dg dg' dg' dg 

dq dp dq dp' 



(25) 
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Eq. ( 22 ) may be rewritten as 

{Aa{q)){t) 
t 

= I dT(^{a{t-T),v,s{q,T)[Af]})^^. (26) 



This is the central result of the paper. The above equa- 
tion has a form similar to the Kubo formula for the re- 
sponse of a dynamical quantity defined in the full 2N- 
dimensional phase space to an external perturbation |13j . 
The relation of formula ( 26 ) with more general ones de- 



rived by Ruelle |16j in the context of dynamical system 
theory remains to be investigated. In the following sec- 
tion, we discuss the special case of a homogeneous QSS, 
i.e., fa{q,p) = P{p) is a function solely of the momentum, 
to obtain an explicit form of the formal solution (21). 



III. HOMOGENEOUS QSS 

We consider a homogeneous QSS with fo{q,p) — P{p), 
where P{p) is any distribution of the momentum, with 
the normalization 



JdqdpPip) = l, JdpP{p) = ^, 



where 



V 



dq 



is the total volume of the coordinate space. 
For a homogeneous QSS, Eq. (12 1 gives 



^Q{q,p)[fo] = -p 



d_ 
dq' 



(27) 



(28) 



(29) 



so that Eq. (21) becomes 
t 

/.fi ,N /■> _(t_^)p^ at;eff(g,T)[A/] dP(p) 
Af(q,p,t) = JdTe "'i^" 



(30) 

which implies that the spatial Fourier and temporal 
Laplace transform of Af{q,p,t) satisfies [17] 

dP{p) 



Afik,p,u;) 



dp 



-ikL[e~""''] 



2TTv{k) J dp' A f{k,p,uj) - K{uj)b{k) 
Here, L denotes the Laplace transform: 

L[e-^P'^]^ dt e'^'^'P'^ ^ ^, 

Jo i(kp-uj) 

assuming Im(a;) to be positive. Thus, 

Afik,p,Lo) 
_ dPjp) k 
dp kp — Lo 



(31) 



(32) 



2TTv{k) I dp' Af{k,p',uj) - K{uj)b{k) 

(33) 



Now, integrating both sides of Eq. (33) with respect to 
p gives 

, -r-,,, , K{uj)b(k) /e{k,uj)~l\ 
Here, e{k,uj) is the "dielectric function" [T]: 

where the integral has to be performed along the Landau 
contour LC that makes Eq. (34 1 valid in the whole of 
the w-plane; we have PT] [T5]: 



l~2nkv{k) J 



e{k,uj) 



kp — LO 

(Im(w) > 0), 



1^2TTkv{k)P J 



-uj dp 



dPjp) 



/k 

(Im(L.) = 0), 



l-2TTkv(k) [ 

^ / J kp—uj op 
— oo 

-iAT,^v{k) ^ 



UJ jk 

(Im(c.) < 0), 

(36) 

where P denotes the principal part. 

From Eq. ( 14 ) , the change in the distribution due to 



the external field is 



/(g,p,t)-P(p) 



1 

2^ 



duj e-"^* J dk e''"^Af{k,p,uj), 

(37) 

where C is the Laplace contour. Integration over p gives 
dp f{q,p,t) = ^ 



2-K 



du e-*"* / dk e''"^ 



K{uj)b{k) (e{k,uj) - 1 



V e{k,uj) ) 



2TTv{k) \ e{k,uj) 



(38) 



where we have used Eq. (34). Let us suppose that the 



expression enclosed by square brackets has singularities 
which are isolated poles of any order. Let {a;p(fc)} be the 
set of poles, while {rp{k)} is the set of residues at these 
poles. Then, by the theorem of residues, we get 



dpf{q,P,t) ^ ^7 + 1^ 

V ZTT 



dk e*'=9^(27ri)rp(A:)e- 



(39) 



From Eq. ([38| , we see that the poles correspond either 
to poles of K{u!) or to the zeros of the dielectric function, 
i.e., values Wp(fc) (complex in general) that satisfy 



e(fc,Wp(fc)) = 0. 



(40) 
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Equation (39) implies that these poles determine the 
growth or decay of the difference J dp f{q,p,t) — y 
in time depending on the location of the poles in the 
complex-w plane. For example, when there are poles in 
the upper-half complex w-plane, the difference grows in 
time. If, on the other hand, the poles lie either on or be- 
low the real-oj axis, the difference does not grow in time, 
but oscillates or decays in time, respectively. 



We have to ensure that our analysis leading to Eq. ( 39 ) 
is consistent with the decomposition in Eq. (|14|) for 



perturbations about a stable stationary state fo{q,p) — 
P{p)- It is thus required that / dp f{q,p,t) — y does 
not grow in time, which means that the aforementioned 
poles cannot lie in the upper-half w-plane. Now, since 
K{t) was chosen to satisfy the conditions K{t = 0) = 
and K{t ^ oo) = a, constant much smaller than 1, it 
follows that K{uj) cannot have poles in the upper-half cj- 
plane. We therefore conclude that Eq. (39) is valid when 



the poles Wp(fc) that come from the zeros of e(k, uj) satisfy 



e{k,ojp{k)) = 0; lm{wp{k))<Q, 



(41) 



corresponding to linear stability of the stationary state 
P{p)- The condition Im(a;p(fc)) = corresponds to 
marginal stability of P{p)- In this case, the zeros of 
the dielectric function lie on the real-w axis so that 
Wp(fc) = Wpr(^) is real. From Eqs. (40) and (36), we 



find that a;pi-(fc) satisfies 



.0 

l-27rfcv(fc)P / 

J —i 



dp dP{p) 
kp — ujpi.{k) dp 



-i2T:^v{k) 



dP{p) 



dp 



= 0. 



(42) 



Equating the real and the imaginary parts to zero, we 
get 



1 - 27rfcw(fc)P 



dp dP{p) 
kp ~ LOpr{k) dp 



= 0, 



dPjp) 
dp 



(43) 
(44) 



Wpr(fe)/fc 



We now move on to apply our analysis to the Hamil- 
tonian mean-field model, a paradigmatic model of long- 
range interactions. 



IV. APPLICATION TO THE HAMILTONIAN 
MEAN-FIELD MODEL 

A. Model 

The Hamiltonian mean-field (HMF) model belongs to 
the class of models with Hamiltonian ([T]), with the addi- 
tional feature that the potential v{q) is even: 



v{q) = 1 ~ cos g. 



(45) 



so that the Hamiltonian is [5] 



N , 

2 



N 



1=1 



S(q,; - q,)] . (46) 



i<j 



The model describes particles moving on a unit circle un- 
der Hamiltonian dynamics ^ . The canonical coordinate 
qi e [0, 27t] specifies the angle for the location of the ith 
particle on the circle with respect to an arbitrary fixed 
axis, while pi is the conjugate momentum 19 . 

The model in the equilibrium state shows a continuous 
transition from a low-energy clustered phase, in which 
the particles are close together on the circle, to a high- 
energy homogeneous phase corresponding to a uniform 
distribution of particles on the circle. The clustering of 
the particles is measured by the magnetization vector 
(m)(t) with components 



i{'m'x){t), {my)it)) ^JJ dqdp {cosq,smq)f{q,p,t), 

(47) 

and mag nitude {m){t) = y/{m.^)'^{t) + {myy{t). In 
terms of {m){t), the energy density is 



(48) 



T of the system: 

r>2 



where the kinetic energy ( ^ ) (0 defines the temperature 



y)W= JJ dqdp^fiq,p,t) 



T 
2' 



(49) 



Note that e is conserved under the dynamics. 

In equilibrium, the single-particle distribution assumes 
the canonical form, fcq{q,p), which is Gaussian in p with 
a non-uniform distribution for q below the transition en- 
ergy density Cc and a uniform one above |20j : 



fcq{q,p) = 



\/l3exp — f3 — TO|*cos(g — 0)^ 



(27r)3/2Jo(/3m? 



(50) 



Here, /q is the modified Bessel function of zero order, f3 
is the inverse temperature, while m^'^ is the equilibrium 
magnetization that decreases continuously from unity at 
zero energy density to zero at Cc and remains zero at 
higher energies. The arbitrary phase (f> in Eq. (50) is 



a result of the rotational invariance of the Hamiltonian 
( 46 ) . The energy at equilibrium is 



1 

2^ 



1 



(51) 



The phase transition in the HMF model occurs within 
both microcanonical and canonical ensembles [HI 1^ . 
Thus, ensemble equivalence, though not guaranteed for 
long-range interacting systems, holds for the HMF model 
jlj. The microcanonical transition energy is Cc = 3/4, 
which corresponds to a transition temperature Tf. = 1/2 
in the canonical ensemble. 
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B. Linear response of homogeneous QSS 

Consider the QSS distribution foiq^p) ~ P{p) whieh 
is homogeneous in coordinate (thus, {nix) fo — ("^y) fa — 
0), but has an arbitrary normahzed distribution for the 
momentum. Here, we study the response of this QSS to 
the external perturbation 



Now, using the fact that (rrix) fo = {my) f„ = 0, 
Eqs. (|60| and ([6T|) imply that 



{mx){t) = 



h 

2^ 



duj 



VjJ V eil.Lv) / 



e(l,w) 



and 



(m,)(i) = 0. 



(65) 



N 



Hc^t = -K{t) ^ cos 



(52) 



which corresponds to the choice 

h[q) = cosq. (53) 
in Eq. ([s]). The specific K{t) we choose is a step function: 
for i < 0, 



It can be proven straightforwardly from the Vlasov equa- 
tion ^ that Eq. ^ holds also in the non-linear re- 
sponse regime {K(t) not necessarily small) for all homo- 
geneous P{p) which are even in p. 

When the zeros of e(l,w) lie only in the lower-half 
complex-oj plane, Eq. (64) gives the time-asymptotic re- 
sponse: 



m = 



/i for t > 0; /K 1. 



(54) 



lim {rrix) (t) 



l-e(l,0) 
e(l,0) 



(66) 



The changes in the magnetization components due to 
the field are 

(Am,)(t) ^ jj dqdp [f{q,p,t) ~ P(p)) cosq 



dp {Af{l,p,u;) + Afi-l,p,Lj)U55) 



This equation implies diverging magnetization at 
e(l,0) — 0, which is clearly not possible as the magneti- 
zation is bounded above by unity. Therefore, in such a 
case, the linear response theory makes an incorrect pre- 
diction. Thus, we rely on formula (66) only when the 
result is much smaller than unity. 



and 



Note that k{ui), given in Eq. (59), has a pole only at 
u! = 0. Following the discussions in Sec. |III[ we thus con- 
clude that the conditions ( 43 ) and ( 44 ) solely determine 



{Amy){t) dqdp (^f{q,p,t) - P{p)^ sing 



2i 



duj e 



c 



the parameters characterizing the distribution P{p) such 
that it is marginally stable. For the HMF model, we 
need to consider only fc = ±1 in these conditions. Since 



dp ( A/(-l,p,w) - A/(l,p, w) )(.56) e(l,w) = e(-l,'^), we write a;pi.(l) = ujp 



that these conditions become 



.(-1) = Wpr, SO 



Using 



5k- 



5k ~i + 5k,i 
2 '■ 



h 



and Eq. (34) in Eqs. (55) and (56 1 gives 
h 



{Amx){t) = i— I duj e 
27r Jc 



V e(l. 



and 



(57) 
(58) 
(59) 

(60) 
(61) 



1 + ttP 



dP{p) 



dp dP{p) 



dp 



-oo P T W] 

= 0. 



pr 



dp 



0, 



(67) 
(68) 



We now consider two representative P{p) and obtain 
the linear response of the corresponding QSS by using 
Eq. (64). For the first case, we obtain the full temporal 



(Am,)(t) = 0. 
Here, we have used the fact that for the HMF model, 

e(l,c^) = e(-l,c^), (62) 



behavior of the response, while in the second case, we 
discuss only the time-asymptotic response. 



1. Water-hag QSS 

The water-bag state corresponds to coordinates uni- 
formly distributed in [0, 27r] and momenta uniformly dis- 
tributed in [—pq,pq]: 



Pip) 



1 1 
27r 2po 



as may be easily checked by using Eq. (57) in Eq. (35) 
It may also be seen that 



0(P + Po) - 0(P-Po) ; J3e[-Po,Po]- 

(69) 

Here, Q{x) denotes the unit step function. The energy 
density is obtained from Eq. (48) as 



e{k,uj) = 1 for fc 7^ ±1. 



(63) 



(70) 
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The dielectric function may be obtained straightfor- 



wardly by using Eq. (351 to get 



e(l,^) 



(71) 



which is analytic in the whole of the w-plane, except at 
the two points uj = ±po- 

As discussed in Sec. |III[ the zeros of the dielectric 
function determine the temporal behavior of the differ- 
ence J dp f{q,p,t) — l/V (here V — 2tt). The zeros of Eq. 

^ occur at Up = ±Vpo " 1/2- For po < p*^ ^ 1/V2, 
(correspondingly, e < e* ~ 7/12), the pair of zeros lies 
on the imaginary-w axis, one in the upper half-plane and 
one in the lower. The one in the upper half-plane makes 
the water-bag state linearly unstable for e < e*. As 
e approaches e* from below, the zeros move along the 
imaginary-w axis and hit the origin when e — e* . At 
higher energies, the zeros start moving on the real-cj axis 
away from the origin in opposite directions. The fact 
that the zeros of the dielectric function are strictly real 
for e > e* implies that the water-bag state is marginally 
stable at these energies, and is therefore a QSS at these 
energies. 

From the discussions in Sec.|III|and those following Eq. 



( 66 1, it follows that the result of the linear Vlasov theory, 
Eq. (64), is valid and phys ically meaningful only when 



Pq > "T/2. Using Eq. ( [flj ) in Eq. (64 1 and performing 
the integral by the residue theorem gives 



{■mx){t) 



2h 



-sin^f^ 
1 V2 



Pi 



pI>\- (72) 



Thus, the linear Vlasov theory predicts that in the pres- 
ence of an external field along x, the corresponding mag- 
netization exhibits oscillations for all times and does not 
approach any time-asymptotic constant value. This pre- 
diction is verified in numerical simulations discussed in 



1 



Sec. VA The average of {mx)(t) over a period of oscil- 
lation is 

1 h 

(ma:)Timc average = ^ J '^^ {f^x) {t) = 2p'^ - V ^°^2 

(73) 

where T is the period of oscillation. In Sec. |Vj we will 
compare this average with numerical results. 



analytic computations of various physical quantities is 
possible. As /3 — cx) the Fermi-Dirac state converges to 
the water-bag state ([69| with po = JJi- 



As shown in the Appendix, to leading order in 1//?^, 
the normalization is given by 



A 



while the energy density is 



1 



24^2^2 



6 



1 



6^2^2 



(75) 



(76) 



satisfies 



dPjp) 
dp 



p=0 



Let us now investigate the conditions (67) and (68) 
for the marginal stability of the state ( 74 ) . Since P(p) 

0, the condition (68) implies that 



0, which on substituting in condition (67) gives 



e(l,0) = 0, (77) 
where, as shown in the Appendix, to order l//3^, we have 

6^V^ 



(78) 



Solving Eq. ( 77 ) gives /i* , the value of /i at the marginal 
stability of the state (74). To order 1//?^, we get 



27r2 
3/32' 



(79) 



which gives the corresponding energy density 



7 

12 



(80) 



such that at higher energies, the state (74) is a QSS 



Following our earlier discussions on the regime of va- 
lidity of the linear Vlasov theory, and using Eq. (78) in 



Eq. (66), we get 



hi 1 



WLx 



6,32^2 



2^-1 



(81) 



C. Linear response of the homogeneous 
equilibrium state 



2. Fermi-Dirac QSS 

We now consider a Fermi-Dirac state in which the co- 
ordinate is uniformly distributed in [0, 27r], while the mo- 
mentum has the usual Fermi-Dirac distribution: 



P{p) = A 



1 



1 



27r 1 + e'3(p'-A') ' 



p G [— oo, oo] 



(74) 



Here, /3 > and /i > are parameters characterizing the 
distribution, while A is the normalization constant. We 



consider the state (74) in the limit of large (3 in which 



It is interesting to consider the response of the distri- 
bution (50) with magnetization = 0, which is the 



equilibrium state of the HMF model for energies e> Cc 
We thus consider the choice 



P{p) 



exp 



2 



(82) 



It is known that the equilibrium state ( 82 ) is also a 
QSS [1 . Indeed, stability condition (68) gives Wpr = 0, 
so that Eq. (67) gives 



e(l,0) = 0, 



(83) 
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where e(l, 0) is given by 



e(l,0) = l-^. 



(84) 



Thus, the state (50) is marginally stable at P — 2, and 
correspondingly, e — e* = 3/4 = Cc- For e > e*, the state 
is a QSS and also the Boltzmann-Gibbs equilibrium state. 



Using Eqs. (84) and ( |66[ ), one gets 

= — ; ; /3 < 2. 

2/^-1' ^ 



(85) 



Therefore, under the perturbation, Eqs. (52) and (54), 



the equilibrium state evolves to an inhomogcneous QSS 
predicted by our linear response theory. Let us compare 
the value of frix in Eq. ( 85 ) with the one predicted by 
equilibrium statistical mechanics, m^'?(/3, h), at the same 
values of the energy and h. This latter quantity is ob- 
tained by solving the implicit equation [T] , 



X 

J 



h{X) 



(86) 



with Ii{X) the modified Bessel function of first order, 
and using the solution X{fi,h) to get 



eq,o u\ ^l(^) 

^q[X) 



(87) 



The corresponding energy is 
2/3 ^ 2 



-hml^{p,h). 



The two values given in Eqs. (85) and (87) are in gen- 
eral different. However, in the high-energy regime, one 
can solve Eq. ( 86 1 for small X to obtain for the equilib- 



rium magnetization the same formula as the one obtained 
by the linear response theory, Eq. ( 85 ) . While compar- 



ing the two magnetization values with numerical results 
at high energies in Sec. |VB[ we are thus not able to dis- 
tinguish between equilibrium and QSS magnetization in 
the presence of the field. 



COMPARISON WITH A^-PARTICLE 
SIMULATIONS 



To verify the analysis presented in Sec. |IV[ we per- 
formed extensive numerical simulations of the A^-particle 
dynamics ^ for the HMF model for large N . The equa- 
tions of motion were integrated using a fourth-order sym- 
plectic scheme [52], with a time step varying from 0.01 
to 0.1. In simulations, we prepare the HMF system at 
time t = in an initial state by sampling independently 
for every particle the coordinate q uniformly in [0, 27r] 
and the momentum according to either the water-bag. 



the Fermi-Dirac, or the Gaussian distribution. Thus, the 
probability distribution of the initial state is 

N 

P{qi,Pi,q2,P2, ■ ■ ■,qN,PN) = YiPiPi) (89) 



where P{p) is given by either Eq. ( 69 ) , ([74h , or ( 82 1 . The 



energy of the initial state is chosen to be such that it is 
a QSS. Then, at time tg > 0, we switch on the external 
perturbation, Eqs. (52) and (54), and follow the time 



evolution of the x-magnetization. 

In obtaining numerical results, two different ap- 
proaches were adopted. In one approach, we followed 
in time the evolution of a single realization of the ini- 
tial state. These simulations are intended to check if 
our predictions based on the Vlasov equation for the 
smooth distribution f(q,p,t) for infinite N are also valid 
for a typical time-evolution trajectory of the empirical 
measure fd{q,p,t) for finite N, where the initial condi- 
tion fd{q,p,0) is obtained from Eqs. (89 1 and ([5]), while 
f{q,p, 0) = P{p)- Rigorous results from Braun and Hepp 
and further analysis by Jain et al. show that these typi- 
cal trajectories stay close to the trajectory of f{q,p, t) for 
times that increase logarithmically with N |9l [15] . When 
P{p) is a stable stationary solution of the Vlasov equa- 
tion, it is known numerically |3] and analytically |23j that 
these times diverge as a power of N , and are therefore 
sufficiently long to allow us to check even for moderate 
values of TV the predictions of our linear Vlasov theory 
for perturbations about Vlasov-stable stationary solution 
P{p). 

In a second approach, we obtained numerical results by 
averaging over an ensemble of realizations of the initial 
state. The time evolution that we get using this second 
method is different from the first one. This approach 
allows us to reach the average and/or asymptotic value 
of an observable, here {mx){t), on a faster time scale 
because of a mechanism of convergence in time, as we 
describe below. 



A. Linear response of homogeneous QSS: Single 
realization 

The oscillatory behavior of {mx){t) predicted for the 
water-bag state, see formula (72), is checked in Fig.[lja). 
Oscillations around a well-defined average persist indef- 
initely with no damping, as predicted by the theory. In 
the inset of the same panel, the theoretical prediction 
is compared with the numerical result for a few oscilla- 
tions. While the agreement is quite good for the first two 
periods of the oscillations, the numerical data display a 
small frequency shift with respect to the theoretical pre- 
diction. Moreover, an amplitude modulation may also 
be observed. We have checked in our A^-particle simula- 
tions that different initial realizations produce different 
frequency shifts, which has a consequence when averag- 
ing over an ensemble of initial realizations, as discussed 
below. 
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In Fig.[IJb), we show {mx){t) for the Fermi-Dirac QSS. 
In this case, we have the theoretical prediction only for 
the asymptotic value rfix given in Eq. (jsTJ . The time evo- 
lution of {mx){t) displays beatings and revivals of oscilla- 
tions around this theoretical value, shown by the dashed 
horizontal line in the figure. There is no sign of asymp- 
totic convergence, even running for longer times. For this 
high value of /3, which makes the Fermi-Dirac distribu- 
tion very close to the water-bag one, we cannot conclude 
that there will be damping in time. We have observed 
a damping for smaller values of /3 when the Fermi-Dirac 
distribution comes closer to a Gaussian. 



B. Linear response of the homogeneous 
equihbrium state: Single realization 

In Fig.[2j we show {mx){t) for the Gaussian QSS. After 
the application of the external field, the magnetization 
sharply increases and fluctuates around a value which is 
slightly below the theoretical prediction, Eq. (85 1. Gon- 



vergence to this latter value is observed on longer times. 



C. Average over initial realizations 

In this section, we present numerical results for the 
three initial QSSs (water-bag, Fermi-Dirac, Gaussian), 
obtained after averaging the time evolution of {nix) (t) 
over a set of realizations (typically a thousand) of the 
initial state. We define the average 



Ensemble average (^-^ 



1 



m,)„(i), (90) 



where (•)„ labels the sample and is the total number 
of different realizations. 

In all cases, we observe a relaxation to an asymptotic 
value. For the water-bag distribution, this value com- 
pares quite well with the time-averaged magnetization 
given in formula (73 1, see Fig. [s] panels (a) and (b). The 
mechanism by which the relaxation to the asymptotic 
value occurs in the water-bag case, in the absence of a 
true relaxation of a single initial realization, is the fre- 
quency shift present in the different initial realizations. 
This leads at a given time to an incoherent superposition 
of the oscillations of the magnetization. For other distri- 
butions, the numerically determined asymptotic value is 
compared with the theoretical value for the single real- 
ization ffix, given in Eq. (81) and Fig. [sjpanels 2(a) and 
2(b), and formula (85) and Fig. [s] panels 3(a) and 3(b). 
The agreement is quite good. 



D. Relaxation of QSS to equilibrium 

For finite values of N, the perturbed HMF system fi- 
nally relaxes to the Boltzmann-Gibbs equilibrium state. 



0.15 



0.075 r. 




200 



0.15 



^ 0.075 r... 



I'l'l'lif.. 



II w 



t^o 



25 



250 

Time t 



500 



FIG. 1. (Color online) {mx){t) vs. time t for the (a) water- 
bag QSS, and (b) Fermi-Dirac QSS, in the HMF model under 
the action of the perturbation, Eqs. (521 and (54 1, with h = 
0.1 switched on at time to ~ 25. (a) The full line in the 
main plot shows the result of A^-particle simulation, while the 
dashed horizontal line is the theoretical time-averaged value 
of {mx){t) given in Eq. (73 1. The system size is A^ = 10^, 
while the parameter po, corresponding to energy e = 0.7, is 
approximately 1.095. In the inset, the numerical result (full 
line) is compared with the theoretical prediction ( 72 1 (dashed 
line), (b) The full line represents simulation results, while 
the horizonta l da shed line is the theoretical asymptotic value 
given in Eq. (81 1. The system size is A = 10^, while /3 = 10 
and n - 



1.2, giving energy 



0.7. 



The presence of a two-step relaxation of the initial water- 
bag QSS with energy e = 0.65, first to the perturbed 
Vlasov state and then to equilibrium, is shown in Fig. |4] 
for increasing system sizes for perturbation, Eqs. ( 52 ) and 
(54), with h = 0.01. The relaxation to the first magneti- 
zation plateau with value ~ 0.125 predicted by the linear 
response theory takes place on a time of 0(1). The final 



10 



0.035 




500 



FIG. 2. (Color online) {mx}{t) vs. time t for the Gaussian 
QSS i n th e H MF m odel under the action of the perturbation, 
Eqs. (521 and (54 1, with h = 0.1 switched on at time to — 25. 
The line made of pluses represents the result of A''-particle 
simulation, while the dashed horizontal line is the theoretical 



asymptotic value given in Eq. (85 1. The system size is 
10'', while P — 0.5, so that the energy e — 1.5. 



relaxation to the equilibrium value of the magnetization 
« 0.42 occurs on a timescale that increases with system 
size, presumably with a power law that remains to be 
investigated further. 



VI. CONCLUDING REMARKS 

In this paper, we studied the response of a Hamilto- 
nian long-range system in a quasistationary state (QSS) 
to an external perturbation. The perturbation couples to 
the canonical coordinates of the individual constituents. 
We pursued our study by analyzing the Vlasov equa- 
tion for the time evolution of the single-particle phase 
space distribution. The QSSs represent stable stationary 
states of the Vlasov equation in the absence of the ex- 
ternal perturbation. We linearized the perturbed Vlasov 
equation about the QSS for weak enough external per- 
turbation to obtain a formal expression for the response 
observed in a single-particle dynamical quantity. For a 
QSS that is homogeneous in the coordinate, we derived a 
closed form expression for the response function. We ap- 
plied this formalism to a paradigmatic model, the Hamil- 
tonian mean-field model, and compared the theoretical 
prediction for three representative QSSs (the water-bag 
QSS, the Fermi-Dirac QSS and the Gaussian QSS) with 
A^-particle simulations for large N. We also showed 
the long-time relaxation of the water-bag QSS to the 
Boltzmann-Gibbs equilibrium state. 
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Appendix A: Normalization, energy density, and 
stability criterion for the Fermi-Dirac distribution 
[Eq. ^] 



Normalization: Consider the distribution Eq. (74 1 
The normalization A satisfies 



dp 



1 -|_ e/3(p^-M) 



1. 



(Al) 



Changing variables and doing an integration by parts, we 
get 



2I3A / dx 



1 -(_ gl3{x-fj.) 



= 1. 



(A2) 



The left hand side may be written in terms of the 
derivative 9/fd (x) / dx of the Fermi-Dirac-like function 
/FD(a;) = 1/[1 + e'^(^-'^)]. We get 



2A 



dx \fx 



X dx 



(A3) 



In the limit of large /3, the derivative dfFu{x)/dx 
approaches the Delta function: lim/j^oo dfFjj{x)/dx — 
—S{x — fi). In this limit, we may expand y/x in a Taylor 
series about /i. 



//i- 



X — fi {x — 



2^ 8^^3/2 



which on substituting in Eq. ( A3 1 gives 
1 



A[2^Io 
Here, 



-Ii- 



1 



(A4) 



= 1. (A5) 



/3-i.oo 



dx 



dy 



dx 

ey 



(A6) 
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FIG. 3. (Color online) Linear response of a water-bag QSS (panels 1(a), 1(b)), a Fermi-Dirac QSS with /B — 10 (panels 2(a), 
2(b)), and the homogeneous equilibrium state (panels 3(a), 3(b)) for the HMF model under the perturbation, Eqs. (52| and 
(54l, with h — 0.01. All simulation data have been averaged over several thousand realizations of the initial state; for details, 
see text. In each case, panel (a) shows the time evolution of the averaged magnetization ("^a:) Ensemble average(^) ^ obtained 
from A'^-particle simula tion s, and its asymptotic approach either to the time average in Eq. ( 73 1 for the water-bag initial state 
or to nix given in Eq. (811 for the Fermi-Dirac QSS, or to nix given in Eq. (851 for the Gaussian QSS. In panel (b), we show 



the A''-particle simulation results for the asymptotic magnetization as a function of energy (the parameter fi in the Fermi-Dirac 
case). The error bars denote the standard deviation of fluctuations around the asymptotic value. The results compare well 
with the theoretical predictions. The system size A'^ is 16, 000 for panels 1(a), 1(b) and panels 2(a), 2(b), and 10, 000 for panels 
3(a), 3(b). 
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which gives 



1? 0.3- 




FIG. 4. (Color online) Two-step relaxation of the water-bag 
QSS toward the Boltzmann-Gibbs equilibrium: {mx){t) vs. 
time t for increasing system size from A'^ = 2000 to N = 64000 
(left to right). Under the perturbation, Eqs. (521 and (541 
with h = 0.01, the water-bag initial QSS with e = 0.65 relaxes 
to an intermediate inhomogeneous QSS with {nix) ~ 0.125 
(lower horizontal dash-dotted line) and then to the equilib- 
rium state with (nix) ~ 0.42 (upper horizontal dashed line). 
The blue thick lines refer to running averages performed to 
smooth out local fluctuations. 



/3 



dx 



{x- fl)(^- 



-dfFujx) 
dx 



ye" 



dy 



ye" 



[i + eyf 



(A7) 



I2 = f3^ dx {x - yiY 
Jq 

- y^e' 



-dfFDjx) 
dx 



13^00 



dy 



y^ey 
(1 + ey)2 



TT 



(A8) 



Thus, to order 1//3 , we find from Eq. (A5| that 

^2 



12/32^3/2 



1, 



A = 



2^V ' 24/32^2 



(AlO) 



Average energy: T he average energy density is ob- 
tained from Eq. ( 48 ) as 



e = A dp 



p2/2 



1 



l + e/J(p^-M) ' 2' ^^^^^ 
Changing variables and doing an integration by parts, we 
get 



dx x^/^ 



<9/fd(2^)\ , 1 



\ dx ) 



3 ./o V dx J ' 2' 



(A12) 



Expanding x^^"^ in a Taylor series about /i and substi- 
tuting in Eq. (A12) give 



A/.3/2 



X2 + .. 



8/32^2 



(A13) 



Using Eq. (AlO), we find that to 0(l//32), the energy 
density is 



e=^(l 



el^^ 6^2 ^2) 



1 



(AM) 



Dielectric function: Using Eqs. (74) and (35 1, we get 

(A15) 



.(1,0) = 1-A/ d 



-^a; V dx 



Expanding l/ \/x in a Taylor series about /i, and sub- 
stituting in Eq. (A15) give 



e(l,0) ^\-A 



Xq II 



3X2 



/Ji 2/3/1^/2 8/32/^5/2 

so that to 0(l//32), we get 

e(l,0) = l-i^(n 



, (A16) 



2/^2 



(A17) 



(A9) where we have used Eq. (AlO). 
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